'''
Created on Nov 23, 2010

@author: oabalbin
'''

import subprocess


# Key: The contigs in the reference index genome (reference.fai), in the dict (reference.dict)
# and dbsnp.vcf should be the same. For this reason no contigs related to the illumina platform should be 
# included if they do not have equivalent in the dbsnp.
# Note that dbsnp132 is already in vcf format. So for any database already in vcf format you cand do

def adjust_dbsnp_format(dbsnp_file):
    '''
    cat ./dbsnp132_00-All.vcf | awk '$0 !~ "_random" {print}' | awk '$0 !~ "_hap" {print}'| awk '$0 !~ "PAR" {print}' | awk '{sub(/MT/, "M", $0); sub(/chr/, "", $0); print}' > dbsnp132_00-All_processed2.vcf
    '''
    tempoutfile=dbsnp_file.replace('.vcf','_processed.vcf')
    
    args=['cat', dbsnp_file,'| awk \'$0 !~ \"_random\" {print}\' | awk \'$0 !~ \"_hap\" {print}\'| awk \'$0 !~ \"PAR\" {print}\' | awk \'{sub(/MT/, "M", $0); sub(/chr/, "", $0); print}\'']
    
    f = open(tempoutfile, "w")
    retcode = subprocess.call(args, stdout=f)
    print retcode
    f.close()
    
    return tempoutfile


# Databases no in vfc format
def build_dbsnp_rod(dirpath,dbsnp_file_name,sort_script_path, genome_index_path,genome_build):
    '''
    It builds a sorted dbsnp file suitabale for using with the gath toolkit 
    # make the hg19 sorted dbSNP record
    # note that we assume header/comment/annotation lines (starting with '#') have been removed from the file already
    cat /seq/references/dbsnp/downloads/snp129_hg18.txt | awk '($2 !~ "_random") && ($2 !~ "_hap")' > tmp.dbsnp.txt
    ~/dev/GenomeAnalysisTK/trunk/perl/sortByRef.pl --k 2 tmp.dbsnp.txt /seq/references/Homo_sapiens_assembly18/v0/Homo_sapiens_assembly18.fasta.fai > dbsnp_129_hg18.rod
    rm tmp.dbsnp.txt
    '''
    thisdbsnp=dirpath+dbsnp_file_name
    tempoutfile= dirpath+'tmp.dbsnp.txt'
    args1=['cat',thisdbsnp,'|','awk','\'($2 !~ "_random") && ($2 !~ "_hap")\'']
    
    f = open(tempoutfile, "w")
    retcode = subprocess.call(args1, stdout=f)
    f.close()
    # Build the rod database
    outfile=dbsnp_file_name.replace('.txt','_'+genome_build+'.rod')
    f = open(outfile, "w")
    args2=[sort_script_path,'sortByRef.pl','--k','2',tempoutfile, genome_index_path]
    retcode = subprocess.call(args2, stdout=f)
    f.close
    retcode = subprocess.call(['rm',tempoutfile])
    
    return outfile
    
#dbsnp_file = '/exds/projects/alignment_indexes/gatk/hg19_illumina/dbsnp132_00-All.vcf'
#adjust_dbsnp_format(dbsnp_file)

def convert_dbsnptohg19(dbsnp_file):
    '''
    Transform the coordinates of vcf file from ncbi reference genome to 
    ucsc. Change chr field from 1 to chr1
    '''
    inputfile = open(dbsnp_file)
    outfile = open(dbsnp_file.replace('.vcf','.processed.vcf'),'w')
    header=True
    for l in inputfile:
        if l.startswith('#') and header:
            fields= l.strip('\n').split('\t')
            if fields[0]=='#CHROM':
                header=False
            outfile.write(l)
            continue
        if not header:
            fields= l.strip('\n').split('\t')
            outfile.write('chr'+l)
            #if fields[1]=='641':
            #    break
    
    inputfile.close()
    outfile.close()


convert_dbsnptohg19('/exds/projects/alignment_indexes/gatk/hg19/dbsnp132_00-All_processed5.vcf')
